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Abstract 

When constructing an algorithm for the numerical integration of a differential equation, 
one should first convert the known ordinary differential equation (ODE) into an ordinary 
difference equation (OAE). Given this difference equation, one can then develop an 
appropriate numerical algorithm. This technical note describes the derivation of two 
such OAEs applicable to a first order ODE. The implicit OAE has the same asymptotic 
expansion as the ODE itself; whereas, the explicit OAE has an asymptotic expansion 
that is similar in structure but different in value when compared with that of the ODE. 

Many physical processes can be represented by systems of first order ODEs of the form 

X Q + P a X a = Q a (for a = 1,2, , N) , (1) 


or equivalently as 

X a = F a [X p ,t ] 

= Q a [Xp,t)- P a {Xp,t\X Q (for a,/3 — 1,2, . . .,N) , (2) 

where the X a are the N independent variables to be solved for, which may be scalar, vector 
or tensor valued in applications. The parameters P a (a scalar), and F a and Q a (of the 
same rank/type as X a ), are functions of the variables Xp and t, in general. If neither 
P<y or Q a depends on Xp for all a,/?, then the system of equations is said to be linear; 
otherwise, it is nonlinear. The dot ‘ ’ ’ is used to denote differentiation with respect to time, 
We choose time to be the dependent variable for illustrative purposes, as it so often is in 
physical applications; however, this is not a necessary restriction on the theory presented 
herein. 

*A technical note prepared for the Journal of Computational Physics. 
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The ODE of equations (1) and (2) is a point function in time. Numerical integration 
algorithms, however, are based on functions evaluated over an interval in time. We therefore 
introduce the first order OAE 



^ = r a [X 0 , At,r), 

(3) 

where 


AX a = X a [t + At]-X a [t], 

(4) 

and therefore 


X a [t + A£] — T F a [Xf3 y At y r]* At . 

(5) 


The OAE represents the ODE in numerical integration algorithms. The formulation is 
implicit when r = t + At; whereas, it is explicit when t = t. The development of existing 
algorithms begins with the tacit assumption that F a = F a . We shall now show that this 
assumption is in error. This is accomplished by analytically integrating the ODE of (1) in 
a recursive manner, which permits the difference function F a in (3) to be determined for 
the OAE. 

One can introduce an integrating factor into the differential equation (1) and thereby 
obtain the recursive integral equation [1] 


rt+At 

: q [<+a<] = exp - / p a [Xi3,z]dz 

h=t 


X a [t\ 


rt+At rt+At 

+ / exp - / Pa[X 0 ,C]d( 

Jt=t J <=( 


Q a [Xp ,£] 


( 6 ) 


which is an exact solution to this first order ODE. John Bernoulli [2] developed a nonre- 
cursive solution similar to (6) in 1697 for the equation X = aX + bX n , whose solution 
was sought by Jacob Bernoulli [3] in 1695. John’s solution was expressed as a quadrature 
because the integral of dzj z in the form of a logarithm was not generally known until later 
that same year [4]. 

The authors have obtained a variety of approximate and exact solutions to (6) by rep- 
resenting the parameters P a and Q a with series expansions, and integrating them term 
by term. The resulting linear approximation, which was acquired by using Taylor series 
expansions, is given by [1] 


X a [f+Ai] 




X 0 [*] + (l 


Q a [ x ^ T ] i n 
P a [X 0 ,r] +U 


d QJX^t] 
dr P a [X 0 ,r] 


• (7) 


This approximation becomes exact whenever P a and Q a are both constants, i.e. whenever 
the ODE is linear. 

Rearranging the recursive solution (7) into a difference equation, one obtains two new 
relationships; they are, for the implicit case, 

(\ - e -PalX fi ,t+At)-At\ 

^ot[Xj3) At , t + AJ] = ( p t + At]- At ) ‘AXa + F a [Xp y t-\- A^]) 

+ of - /o\ 

d(t-\-At) P a [Xp, t + At] 
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and for the explicit case, 


F a [Xp, At , t] 


l _ e -PalX fl ,t]-±t\ 

Pa[X p,t]' At J 


F a [Xp y t] + 0 


d Q«[Xp,t] 

dtP a [Xp,t) 


(9) 


Notice the presence of an extra AX a term in the implicit relation (8), which acts as a 
correction to the derivative F Q , and which is not present in the explicit relation (9). 

The coefficient (1 — e~ p< *[ x P' T ]' At )/ P a [Xp> r]-At is a correction factor of the first order 
that results from taking the differential function F a at time r and converting it into a 
difference function T Q at time r taken over the interval (t,t+At). (Higher order solutions 
for the difference function will be given in a future paper.) This coefficient goes to 1 in 
the limit as At goes to 0, as it must so as to recover the differential; in other words, 


Yim^ Q \X 0 ,At,T] = F a [X 0 , t] 


( 10 ) 


for both the implicit and explicit cases. A note of caution when writing computer code. In 
the neighborhood of P Q [Xp, r]-At « 0, one needs to expand (1 — e~ p<3, ^ x ^* T ^^ i )/P a [Xp 7 r]-At 
into a power series to secure a sound computational algorithm. 

The presence of the coefficient (1 P a [Xp, r]-At also introduces desirable 

asymptotic characteristics into our numerical approximations. In particular, for the implicit 
case where r = t + A t, the asymptotic expansion for the OAE (3 with 8) is given by 


lim X a [t + At] 

At — ►large 
Pa> 0 


Qg[Xp,t + At] 
P a [Xp,t + At\ ’ 


( 11 ) 


which is also the asymptotic expansion of the ODE (1 & 2). In contrast, for the explicit 
case where r = t, the asymptotic expansion for the OAE (3 with 9) is given by 


lim X a [t + At] 

At— + large 
Pa>0 


QcXxp, t] 
P a [X 0 ,t) ■ 


(12) 


These two asymptotic expansions differ only in when their parameters P a and Q a are 
evaluated. For large time steps, the implicit case is asymptotically accurate and stable for 
exponentially decaying solutions, i.e. when P a [Xp,t\ > 0 Vi. Stability becomes an issue 
only when P a < 0. In constrast, the explicit case will oscillate around the true solution for 
large time steps, but with much less potential of becoming unbounded when compared with 
equivalent algorithms constructed without our correction coefficient. These oscillations can 
be mitigated only by choosing smaller time steps. 

The second integral in (6) is a Laplace [5] integral where the integrand has its largest 
value at the upper limit, t-\- A£, and possesses an evanescent memory of the forcing function, 
Q a [X0,f], provided that P a [Xp,(] > 0 over the interval (/,t + A/). This fading memory 
means that the solution will depend mainly on the recent values of the forcing function, 
and that by concentrating the accuracy on the recent past we obtain accurate asymptotic 
representations of the solution. 

In the implicit solution, the integrands in (6) were expanded in Taylor series about their 
upper limits [1], where each integrand has its largest value and contributes the most to the 
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integral. By retaining but a single term in the Taylor series expansions the integrands are 
accurately approximated where they are largest, and the neglect of the higher order terms 
is only felt near the lower limits where each integrand contributes only a small amount to 
the integral because of its exponential decay from the upper limit. The neglect of the higher 
order terms in the Taylor series thus results in an algorithm that is asymptotically correct 
at the upper limit. Normally, when treating asymptotic expansions, the exponential decay 
of the integrand allows the lower limit to be replaced with zero or minus infinity to ease the 
integration. This was not done in the present case, however, so that by retaining the lower 
limit as t, we obtain a uniformly valid asymptotic algorithm in the implicit approximation 
provided that P a [X@,(] > 0 over the interval (t,t + At). 

In the explicit solution, the integrands in (6) were expanded in Taylor series about their 
lower limits [1], where the neglect of the higher order terms in the Taylor series results in 
integrands that become progressively more inaccurate as they approach their upper limits 
where the contribution from each integrand is most important. The explicit approximation 
is not, therefore, a valid asymptotic representation of the integral when the Taylor series 
is truncated at a finite number of terms. However, when P a [Xp,£] < 0 over the interval 
( t,t+At ), the reverse situation occurs. In this case the asymptotic solution is now obtained 
by expanding the integrands about their lower limits, where this region of each integrand 
now contributes the most to the integral. Here the implicit method, obtained by expanding 
about the upper limit, does not give a valid asymptotic representation of the integral. 

In conclusion, one only needs to replace the differential function F a [Xp, r] with the ap- 
propriate difference function T a [Xp, At, r] in many existing numerical integration methods 
( e.g . Euler and Runge-Kutta) to construct an appropriate OAE for the numerical integra- 
tion of a given ODE, and thereby obtain substantial improvements in their performance. 
We have demonstrated this in references [1, 6]. 
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